■■■■ 


■ 


■ 


■ 


£9 


■ 


■ 


4*2 


■ 


■l 


■ 


■ 


■ 


■ 


■ 


SUB 
■  ■■ 


;•>' 


MITED  STATES 
'AflTMENT  OF 

IMMERCE 
IBLICATION 


* 


NBS  TECHNICAL  NOTE  625 


*>•««•  0« 


U.S. 

ARTMENT 
OF 
IDMMERCE 

National 
Bureau 
of 
Standards 


Computer  Programs  for  Thermodynamic 
and  Transport  Properties  of  Hydrogen 

(Tabcode-ll) 


NATIONAL  BUREAU  OF  STANDARDS 

The  National  Bureau  of  Standards1  was  established  by  an  act  of  Congress  March  3, 
1901.  The  Bureau's  overall  goal  is  to  strengthen  and  advance  the  Nation's  science  and 
technology  and  facilitate  their  effective  application  for  public  benefit.  To  this  end,  the 
Bureau  conducts  research  and  provides:  (1)  a  basis  for  the  Nation's  physical  measure- 
ment system,  (2)  scientific  and  technological  services  for  industry  and  government,  (3) 
a  technical  basis  for  equity  in  trade,  and  (4)  technical  services  to  promote  public  safety. 
The  Bureau  consists  of  the  Institute  for  Basic  Standards,  the  Institute  for  Materials 
Research,  the  Institute  for  Applied  Technology,  the  Center  for  Computer  Sciences  and 
Technology,  and  the  Office  for  Information  Programs. 

THE  INSTITUTE  FOR  BASIC  STANDARDS  provides  the  central  basis  within  the 
United  States  of  a  complete  and  consistent  system  of  physical  measurement;  coordinates 
that  system  with  measurement  systems  of  other  nations;  and  furnishes  essential  services 
leading  to  accurate  and  uniform  physical  measurements  throughout  the  Nation's  scien- 
tific community,  industry,  and  commerce.  The  Institute  consists  of  a  Center  for  Radia- 
tion Research,  an  Office  of  Measurement  Services  and  the  following  divisions: 

Applied  Mathematics — Electricity — Heat — Mechanics — Optical  Physics — Linac 
Radiation2 — Nuclear  Radiation2 — Applied  Radiation2 — Quantum  Electronics3 — 
Electromagnetics3 — Time  and  Frequency3 — Laboratory  Astrophysics3 — Cryo- 
genics3. 

THE  INSTITUTE  FOR  MATERIALS  RESEARCH  conducts  materials  research  lead- 
ing to  improved  methods  of  measurement,  standards,  and  data  on  the  properties  of 
well-characterized  materials  needed  by  industry,  commerce,  educational  institutions,  and 
Government;  provides  advisory  and  research  services  to  other  Government  agencies; 
and  develops,  produces,  and  distributes  standard  reference  materials.  The  Institute  con- 
sists of  the  Office  of  Standard  Reference  Materials  and  the  following  divisions: 

Analytical     Chemistry — Polymers — Metallurgy — Inorganic     Materials — Reactor 

Radiation — Physical  Chemistry. 

THE  INSTITUTE  FOR  APPLIED  TECHNOLOGY  provides  technical  services  to  pro- 
mote the  use  of  available  technology  and  to  facilitate  technological  innovation  in  indus- 
try and  Government;  cooperates  with  public  and  private  organizations  leading  to  the 
development  of  technological  standards  (including  mandatory  safety  standards),  codes 
and  methods  of  test;  and  provides  technical  advice  and  services  to  Government  agencies 
upon  request.  The  Institute  also  monitors  NBS  engineering  standards  activities  and 
provides  liaison  between  NBS  and  national  and  international  engineering  standards 
bodies.  The  Institute  consists  of  the  following  divisions  and  offices: 

Engineering  Standards  Services — Weights  and  Measures — Invention  and 
Innovation — Product  Evaluation  Technology — Building  Research — Electronic 
Technology — Technical  Analysis — Measurement  Engineering — Office  of  Fire 
Programs. 

THE  CENTER  FOR  COMPUTER  SCIENCES  AND  TECHNOLOGY  conducts  re- 
search and  provides  technical  services  designed  to  aid  Government  agencies  in  improv- 
ing cost  effectiveness  in  the  conduct  of  their  programs  through  the  selection,  acquisition, 
and  effective  utilization  of  automatic  data  processing  equipment;  and  serves  as  the  prin- 
cipal focus  within  the  executive  branch  for  the  development  of  Federal  standards  for 
automatic  data  processing  equipment,  techniques,  and  computer  languages.  The  Center 
consists  of  the  following  offices  and  divisions: 

Information  Processing  Standards — Computer  Information — Computer  Services 
— Systems  Development — Information  Processing  Technology. 

THE  OFFICE  FOR  INFORMATION  PROGRAMS  promotes  optimum  dissemination 
and  accessibility  of  scientific  information  generated  within  NBS  and  other  agencies  of 
the  Federal  Government;  promotes  the  development  of  the  National  Standard  Reference 
Data  System  and  a  system  of  information  analysis  centers  dealing  with  the  broader 
aspects  of  the  National  Measurement  System;  provides  appropriate  services  to  ensure 
that  the  NBS  staff  has  optimum  accessibility  to  the  scientific  information  of  the  world, 
and  directs  the  public  information  activities  of  the  Bureau.  The  Office  consists  of  the 
following  organizational  units: 

Office  of  Standard  Reference  Data — Office  of  Technical  Information  and 
Publications — Library — Office    of    International    Relations. 


1  Headquarters  and  Laboratories  at  Gaithersburg,  Maryland,  unless  otherwise  noted;  mailing  address  Washing- 
ton, D.C.  20234. 

2  Part  of  the  Center  for  Radiation  Research. 

3  Located  at  Boulder,  Colorado  80302. 


AL  BU21E&U  OF  STANDARDS 

jftM  2  9  1S73 

lc"°°    Computer  Programs  for  Thermodynamic  and 


: 


53 

Transport  Properties  of  Hydrogen  (Tabcode-ll) 


H.  M.  Roder 
R.  D.  McCarty 
W.  J.  Hall 


Cryogenics  Division 
Institute  for  Basic  Standards 
0  i  <  National  Bureau  of  Standards 
Boulder,  Colorado    80302 


NBS  Technical  Notes  are  designed  to  supplement  the 
Bureau's  regular  publications  program.  They  provide 
a  means  for  making  available  scientific  data  that  are 
of  transient  or  limited  interest.  Technical  Notes  may 
be  listed  or  referred  to  in  the  open  literature. 


~&i 


9<"»CAU    0* 


S? 


U.S.   DEPARTMENT  OF   COMMERCE,   Peter  G.  Peterson,  Secretary 
NATIONAL  BUREAU  OF  STANDARDS,   Lawrence  M.  Kushner,  Acting  Director 

Issued  October  1972 


National  Bureau  of  Standards  Technical  Note  625 

Nat  Bur.  Stand.  (U.S.),  Tech.  Note  625,  226  pages  (October  1972) 
CODEN:    NBTNAE 


For  sale  by  the  Superintendent  of  Documents,  U.S.  Government  Printing  Office,  Washington,  D.C.   20402 

(Order  by  SD  Catalog  No.  C13.46:625). 


CONTENTS 

Page 

1 .  Introduction 1 

1.  1       Symbols  and  Units  Used 2 

2.  Hydrogen  Modifications,   Property  Differences,  and  Phase  Diagrams 4 

2.  1       Hydrogen  Modifications 4 

2.2  Property  Differences      5 

2.3  Range  of  the  Programs  and  Phase  Diagrams 6 

2.4  The  Phase  Designator,  Q 8 

3.  Overview  of  the  Programs 9 

3.  1       Method  Used 10 

3.2  Peculiarities  of  the  P-T  Programs 10 

3.3  Peculiarities  of  the  P-H  Programs 10 

3.4  The  Saturation  Routines 11 

4.  Programming  Information 12 

4.  1       Deck  Description 12 

4.2  Calling  Sequence 15 

4.3  Exceeding  the  Range  of  the  Programs 16 

4.4  Conversion  to  Other  Systems  or  Installations 16 

4.5  Pitfalls,   Common  Errors  in  Applications 17 

4.6  Debugging  Aids 17 

5.  Structural  Details  of  the  Programs 18 

5.  1       Grid  and  Region  Layout 18 

5.2  Speed 20 

5.3  Continuous  Functions 20 

5.4  Modifications,   Future  Developments 21 

6.  Estimation  of  Errors 22 

7.  Sources  of  Properties  Data      24 

7.  1       Ideal  Gas  Properties 24 

7.2  The  PVT  Surface  at  Low  Temperatures 26 

7.3  The  PVT  Surface  at  Intermediate  Temperatures 26 

7.4  The  PVT  Extrapolation      27 

7.5  Dissociation 29 

7.6  Viscosity 30 

7.7  Thermal  Conductivity 31 

7. 8  Generating  Functions       33 

8.  References 34 


Page 

Appendix  A.        Tables  of  Values  for  Selected  Isobars 37 

Appendix  B.        Program  Listings:     CDC  and  IBM 52 

Appendix  C.        Plots  of  Maximum  Interpolation  Error  and  Estimated 

Source  Error 177 

LIST  OF  TABLES 

1.  Values  of  the  Phase  Designator,  Q 8 

2.  List  of  Available  Program  Decks 9 

3.  Returns  from  the  Liquid- Vapor  Two  Phase  Region      12 

4.  Deck  Description,   p-hydrogen 13 

5.  Deck  Description,    e-hydrogen 14 

6.  Low  Temperature  Error  Grid,   Pressures,  psia       23 

7.  Low  Temperature  Error  Grid,   Temperature,   °R 24 

8.  Coefficients  of  the  Equation  of  State 27 

9.  Generating  Functions 33 

LIST  OF  FIGURES 

1.  Hydrogen  Composition  at  Equilibrium 4 

2.  Ideal  Gas  Specific  Heat  of  Para  and  Equilibrium  Hydrogen 6 

3.  The  P-T  Phase  Diagram  for  Hydrogen 7 

4.  The  P-H  Phase  Diagram  for  Parahydrogen      7 

5.  Region  Layout  for  C     with  Entry  Variables  P-H,   p-hydrogen 18 

P 

6.  Partial  Grid  Layout  for  Program  PTENTH 19 

7.  Matching  of  Values  at  Boundaries      21 

8.  Primary  Sources  of  Data 25 

9.  Comparison  of  PVT  and  Derived  Properties 28 

10.  The  Effect  of  Dissociation  on  Enthalpy 30 

11.  The  Viscosity  of  Dissociating  Hydrogen  at  Various  Pressures 31 

12.  The  Thermal  Conductivity  of  Dissociating  Hydrogen  at  Various  Pressures  ....  32 


IV 


COMPUTER  PROGRAMS  FOR  THERMODYNAMIC 

AND  TRANSPORT  PROPERTIES  OF  HYDROGEN 

(Tabcode-II) 

H.   M.   Roder,  R.    D.    McCarty,  and  W.   J.   Hall 

The  thermodynamic  and  transport  properties  of  para  and 
equilibrium  hydrogen  have  been  programmed  into  a  series  of  com- 
puter routines.     Input  variables  are  the  pair's  pressure-temperature 
and  pressure- enthalpy.     The  programs  cover  the  range  from  1  to 
5000  psia  (34  MN/m^)  with  temperatures  from  the  triple  point  to 
6000°R  (3300  K)  or  enthalpies  from  -130  BTU/lb  (-623  J/mol)  to 
25,000  BTU/lb  (117000  J/mol).     Output  variables  are  enthalpy  or 
temperature,   density,   entropy,   thermal  conductivity,  viscosity, 
velocity  of  sound,  heat  capacity  at  constant  pressure,   heat  capacity 
at  constant  volume,   the  heat  capacity  ratio,  and  a  heat  transfer 
parameter.     Property  values  on  the  liquid  and  vapor  boundaries  are 
conveniently  obtained  through  two  small  routines.     The  programs 
achieve  high  speed  by  using  linear  interpolation  in  a  grid  of  precom- 
puted  points  which  define  the  surface  of  the  property  returned.     The 
maximum  errors  arising  from  the  linear  interpolation  are  shown  on 
individual  deviation  plots  for  each  combination  of  variables.     Error 
estimates  for  the  sources  of  data  are  similarly  displayed. 

Key  Words:     Computer  programs;  density;  enthalpy;  entropy; 
heat  capacity  at  constant  pressure;  heat  capacity  at  constant 
volume;  heat  capacity  ratio;  heat  transfer  coefficient;  hydrogen; 
pressure;  saturation  boundary;  temperature;  thermal  conductivity; 
velocity  of  sound;  viscosity. 

1.     INTRODUCTION 

Extensive  use  of  hydrogen  in  the  chemical  industry  and  as  a  propellant  in  the  U.    S. 

space  program  requires  that  the  engineer  be  provided  values  of  thermodynamic  and  transport 

properties  of  this  fluid  in  useful  form.     In  many  problems  the  properties  have  to  be  evaluated 

at  numerous  thermodynamic  state  points.     The  problems  are  often  handled  on  a  computer 

and  they  involve  iterative  computational  methods.     In  the  present  report  we  describe  a 

series  of  computer  routines  which  are  fast,   moderately  accurate,  and  operate  over  wide 

ranges  of  input  variables.     They  are  based  on  selected  source  data.     The  input  (entry)  and 

output  (returned)  variables  are: 

Input  Output  Input  Output 

pressure-  enthalpy  pressure-         temperature,  and  phase  or  state 

temperature      density  enthalpy      density 

entropy  entropy 

thermal  conductivity  thermal  conductivity 

viscosity  viscosity 

velocity  of  sound  velocity  of  sound 

heat  capacity,     C  heat  capacity,     C 

heat  capacity,     C  heat  capacity,     C 

heat  capacity  ratio  heat  capacity  ratio 

heat  transfer  parameter  heat  transfer  parameter 


In  addition  there  are  two  small  saturation  routines  which  return  values  on  the  satura- 
tion boundaries.  The  computer  routines  were  developed  primarily  for  nuclear  rocket  engine 
design.     The  programs     cover  a  wide  range  of  conditions     as  follows:    for  temperatures 

from  the  triple  point  (24.845°R,    13.803  K)  to  about  6000°R  (3300  K),   for  pressures  from  1  to 

2 
5000  psia  (34  MN/m   ),  and  for  enthalpies  from  -130  to  25,000  BTU/lb  (-623  to  117000  J/mol). 

In  addition,   the  thermodynamic  and  transport  properties  of  both  para  hydrogen  and  equilibrium 

hydrogen  are  available  through  these  programs.     An  important  criterion  for  many  applications 

is  computational  speed,  as  some  process  control  calculations  have  to  be  made  in  real  time. 

The  linear  interpolation  employed  achieves  speed  but  sacrifices  accuracy.     As  an  aid  to  the 

user,  we  have  determined  the  maximum  error  due  to  the  interpolation  scheme  and  estimated 

the  error  in  the  source  data.     These  errors  are  displayed  in  departure  graphs  for  each 

program. 

An  older  version  of  the  para  hydrogen  programs,   which  has  been  distributed  widely, 

should  now  be  considered  obsolete  (Tabcode-I,   Hall,   et  al.  ,    1967,  and  McCarty  1968).     The 

changes  and  improvements  over  the  older  version  are  as  follows: 

1.  The  temperature  range  has  been  extended  to  6000° R,   the  enthalpy  range  to  25,000 
BTU/lb. 

2.  New  measurements  of  thermal  conductivity  between  30°   and  275°R  are  now  included. 

3.  Improved  values  for  both  thermal  conductivity  and  viscosity  have  been  inserted  for 
temperatures  above  1200°R. 

4.  The  saturation  subroutines  are  now  included  in  the  package. 

5.  For  equilibrium  hydrogen,   which  was  not  treated  previously,   an  entire  set  of  programs 
is  now  provided. 

6.  A  combination  of  variables  is  now  included  which  describe  a  commonly  used  heat 
transfer  parameter,   designated  herein  as  L-factor  (see  Section  1.1). 

The  main  divisions  of  the  report  are  text,  program  listings,  and  error  graphs.     The 
text  contains  all  the  normal  description  given  for  computer  programs  such  as  deck  size, 
calling  sequence,    etc.     However,  for  this  problem  the  selection  of  source  data  and  the  com- 
putation of  the  property  values  placed  in  data  arrays  are  important;  they  are,   therefore, 
described  in  some  detail. 

1.  1     Symbols  and  Units  Used 

The  programs  were  developed  for  engineering  applications;  therefore,   the  units  of 
this  report  are  those  commonly  used  in  engineering,    i.e.,  psia,   BTU,    etc.     The  sources 
of  data,  however,   normally  employ  some  version  of  metric  units.     In  the  discussion  of  the 
equation  of  state,   the  virial  B,   and  in  figures  11  and  12,   the  units  of  the  original  sources 
had  to  be  retained.     Conversion  factors  for  changes  from  source  papers  to  generating 
functions  to  the  present  program  have  all  been  taken  from  NASA  Special  Publication  7012 
(Mechtly  196  9). 
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P  pressure,    psia 

T  temperature,    °R 

H  enthalpy,    BTU/lb 

H  enthalpy  of  saturated  vapor,    BTU/lb 

gas 

H,.       .  ,      enthalpy  of  saturated  liquid,    BTU/lb 
liquid 

S  entropy,    BTU/lb  -   °R 

p  density,    lb /ft 

C  specific  heat  at  constant  pressure,    BTU/lb  -   °R 

P 

C  specific  heat  at  constant  volume,    BTU/lb  -  °R 

v 

Y  the  ratio  C    /C    ,    dimensionless 

P      v 

uo  sound  velocity,    ft/s 

T\  viscosity,    lb-h/fts 

k  thermal  conductivity  BTU/h-ft  -  °R 

L  a  heat  transfer  parameter  equal  to  k  "     C°*4/ri  "4,    with  units  as 

P 
shown  above. [  This  heat  transfer  parameter  is  the  combination  of 

thermal  properties  occurring  in  the  conventional  Dittus-Boelter 

(Nusselt  type)  equation  describing  the  forced  convection  heat 

transfer  coefficient  for  turbulent  flow  (see,  for  example, 

equation  9- 10a,  page  219,   McAdams  1954)]. 

Q  the  phase  designator,    dimensionless 

Additional  symbols  used  in  text: 

B  or  B(T)  the  second  virial  coefficient 

E°  ideal  gas  internal  energy  at  absolute  zero,    the  reference  state, 

taken  to  be  zero 

tt  general  symbol  for  property 

P  critical  pressure 

c 

R  gas  constant 

RE  general  symbol  for  the  return  from  the  liquid-vapor  two  phase  region 

superscript  °  ideal  gas  thermodynamic  properties,    i.  e.  ,    C°,    C°,    H° 

U  internal  energy 

V  volume 
Z  PV/RT 

Important  state  points  for  para  and  equilibrium  hydrogen: 

triple  point  normal  boiling  point  critical  point 

P  =  1.  022  psia  P  =   14.  696  psia  P  =   187.  506  psia 

T  =  24.  845°R  T  =   36.  482°R  T  =  59.  356°R 


2.     HYDROGEN  MODIFICATIONS,   PROPERTY  DIFFERENCES,  AND  PHASE  DIAGRAMS 

A  brief  review  of  the  more  important  features  in  the  physics  of  hydrogen  is  essential 
if  we  are  to  understand  what  properties  the  programs  return,  and  why  some  of  the 
peculiarities  of  the  programs  occur. 

2.  1    Hydrogen  Modifications 

Hydrogen  molecules  occur  in  two  natural  modifications,  parahydrogen  and  ortho- 
hydrogen.     In  the  parahydrogen  molecule  the  two  nuclear  spins  are  opposed  while  in  ortho- 
hydrogen  they  are  alined.     Under  normal  conditions,   the  transition  probability  between 
states  is  practically  zero;  a  catalyst,  a  sharp  gradient  in  magnetic  field,   or  external  radia- 
tion (e.g.,  nuclear  radiation,  x-rays,  ultraviolet,   etc.)  is  required  to  let  the  transition  pro- 
ceed.    Conversion  between  these  two  modifications  is  temperature  dependent,  as  shown  in 
figure  1,   which  is  based  on  Woolley,   et  al.  ,    (1948).     At  low  temperatures  in  the  presence  of 
an  appropriate  catalyst,   the  thermal  equilibrium  favors  the  para  modification.     At  the 
normal  boiling  point  (36.  482° R),   the  equilibrium  composition  is  99.79%  para  and  only  0.21% 
orthohydrogen.     This  composition  is  often  simply  referred  to  as  parahydrogen  or  p-hydrogen. 


100 


z 

LlJ 

ID 
O 
(Z 

>    75 

X 

< 
< 

Q_ 
UJ 

< 

I- 

U    50 

O 

QL 

LU 

Q_ 


25 


TEMPERATURE,    K 

100  200 


300 


100 


200  300 

TEMPERATURE, °R 


400 


500 


Figure  1.     Hydrogen  Composition  at  Equilibrium 
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Conversely,    the  equilibrium  composition  near  ambient  temperature  is  75%  ortho  and  25% 
para.     This  latter  composition  is  usually  designated  normal  or  n-hydrogen.     When  conversion 
occurs  between  the  normal  boiling  point  and  room  temperature,    then  the  thermodynamic 
equilibrium  mixture  of  para  and  ortho  hydrogen  is  called  "equilibrium"  or  e-hydrogen  and 
the  composition  at  intermediate  temperatures  varies  as  shown  in  figure  1.     In  this  report  we 
describe  programs  for  both  p-hydrogen  and  e-hydrogen. 

At  temperatures  above  2700°R   (1500  K)  the  hydrogen  molecule  dissociates  into  hydro- 
gen atoms.     This  dissociation  is  both  temperature  and  pressure  dependent.     Assumptions 
made  in  this  report  are  that  dissociation  occurs  according  to  the  equilibrium  constant  for 
normal  hydrogen,  and  that  the  property  values  can  be  obtained  by  adding  the  dissociation 
effects  calculated  for  the  ideal  gas  to  the  extrapolated  properties  of  the  real  gas. 

2.2    Property  Differences 

At  comparable  state  points,  property  values  of  the  various  modifications  of  hydrogen 
may  differ,    that  is,    some  of  the  properties  depend  on  ortho-para  composition.     Experi- 
mentally these  effects  are  observed  in  different  vapor  pressures  and  in  different  critical 
parameters.     Most  important  to  this  discussion  is  the  fact  that  the  properties  of  the  ideal 
gas  for  para  and  equilibrium  hydrogen  as  calculated  from  statistical  mechanics  differ,    and 
further,    since  the  composition  change  depends  only  on  temperature,   that  the  property 
differences  depend  only  on  temperature.     Thus,    if  a  difference  exists  for  an  ideal  gas  prop- 
erty at  low  pressures,   then  the  same  difference  exists  at  all  other  pressure  levels  for  a 
given  temperature.      Stated  another  way,    the  statistical  calculation,   by  including  the  temper- 
ature dependent  composition  change,   automatically  includes  the  reaction  occurring  as  para- 
transforms  to  equilibrium- hydrogen  —   i.  e.  ,    the  sensible  heat  of  e-hydrogen  includes  the 
appropriate  part  of  the  heat  of  chemical  reaction.     Not  all  properties  are  composition  depen- 
dent; density  and  viscosity,   for  example,   are  nearly  independent  of  composition.      On  the 
other  hand,   the  properties  with  significant  ortho- para  dependency  are  specific  heat,   and 
properties  related  to  specific  heat  such  as  velocity  of  sound,    entropy,    enthalpy,   and  thermal 
conductivity.     When  the  properties  of  p-  and  e-hydrogen  differ,   they  may  differ  considerably. 
For  example,   the  value  of  enthalpy  near  the  critical  point  for  p-hydrogen  is  about  16  BTU/lb, 
whereas  the  value  for  equilibrium  hydrogen  is  about  31  BTU/lb.     The  difference  in  prop- 
erties for  the  ideal  gas  specific  heats,     C    °  ,     is  illustrated  in  figure  2.     Differences  for  the 

p 

other  properties  are  best  expressed  in  terms  of  a  table;  approximate  differences  can  be 
obtained  by  comparing  corresponding  entries  for  a  pressure  of  1   psia  in  Appendix  A  (pages 
38  and  45). 
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Figure  2.     Ideal  Gas  Specific  Heat  of  Para  and  Equilibrium  Hydrogen 


2.3    Range  of  the  Programs  and  Phase  Diagrams 

For  programs  which  use  pressure  and  temperature  as  entry  variables,   the  pressure 
range  is  1  to  5000  psia  and  the  temperature  range  is  from  the  triple  point  (24.845°R)  to 
6000°R.     The  P-T  range  of  the  programs  for  both  para  and  equilibrium  hydrogen  is  illus- 
trated by  the  shaded  area  in  figure  3.     The  programs  cover  property  values  for  liquid,   gas  or 
vapor  states,   and  dissociation.     Note  that  the  programs  do  not  cover  the  solid  phase  nor 
values  at  pressures  less  than  1  psia.     The  low  pressure  limit  is  arbitrary;  it  arises  from  the 
behavior  of  entropy  as  follows.     At  any  given  temperature,   the  value  of  entropy  approaches 
infinity  as  the  pressure  approaches  zero.     Since  we  cannot  represent  an  infinity  with  linear 
interpolation  between  fixed  points,  a  cutoff,    1  psia,   was  chosen. 

For  the  programs  which  use  pressure  and  enthalpy  as  entry  variables,   the  pressure 
range  is  the  same,    1  psia  to  5000  psia.     The  range  of  enthalpy  is  -130  BTU/lb  to  25,000 
BTU/lb.     The  enthalpy  limits  correspond  to  the  enthalpy  of  the  liquid  at  the  triple  point,   and 
the  ideal  gas  enthalpy  at  the  highest  temperature,    6000°R.     The  phase  diagram  of  para- 
hydrogen  in  pressure- enthalpy  coordinates  is  shown  in  figure  4.     The  phase  diagram  for 
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Figure  3.     The  P-T  Phase  Diagram  for  Hydrogen 
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e-hydrogen  would  look  quite  similar  to  figure  4  except  that  the  liquid- vapor  two  phase  bound- 
ary is  gradually  displaced  to  higher  values  of  enthalpy.     For  example,   the  enthalpy  of 
e-hydrogen  at  the  critical  pressure  is  about  31  BTU/lb  while  the  maximum  value  of  enthalpy 
along  the  two-phase  boundary  is  about  110  BTU/lb.     The  programs  cover  property  values  for 
liquid,   gas  or  vapor  states,   and  dissociation  as  before.     They  do  not  cover  the  solid  phase, 
pressures  below  1  psia,   or  the  two  phase  liquid- solid  ("slush")  conditions. 

A  comparison  of  figures  3  and  4  discloses  that  the  conditions,   the  phases  or  states 
of  the  hydrogen,    covered  by  the  P-T  and  the  P-H  programs  are  similar,    except  that  the  P-H 
programs  cover  the  range  of  states  where  liquid  and  gas  coexist  while  the  P-T  programs  do 
not.     A  concise  designation  of  the  phase  or  state  of  any  given  point  is  therefore  available 
only  by  reference  to  P-H  coordinates.     Since  the  description  of  phase  is  important  in  many 
problems,   a  scheme  to  designate  the  phase  was  devised  and  is  presented  in  the  next  section. 


2.4    The  Phase  Designator,   Q 

The  various  physical  states  possible  in  the  P-H  phase  diagram  are  described  by  the 
phase  designator,  Q,   available  only  from  program  PHTEMP  (see  also  page  11).     Values  of 
Q    are  shown  in  figure  4  and  in  table  1.     Note  that  figure  4  includes  scale  changes  at 
100  BTU/lb  and  at  200  psia.     For  the  special  case  of  the  two  phase  liquid-vapor  region,     Q 
corresponds  to  the  conventional  quality,   as  defined  on  page  11. 


Table  1.     Values  of  the  Phase  Designator,  Q 


Value  of  Q 

Applicable  portions  of  the  phase  diagram,     P       =     187.506  psia 

c 

1.0 

single  phase  gas  or  vapor 

single  phase  compressed  fluid  at  pressures  above  P 

c 
saturated  vapor  line 

0.0  <Q  <  1.0 

two  phase  liquid-vapor  region 

0.0 

saturated  liquid  line 

compressed  liquid  (fluid)  states  at  pressures  below  P 

liquid- solid  boundary  line 

-1.0  <  Q  <  0.0 

two  phase  liquid- solid  region 

-1.0 

solid  -   compressed  solid  boundary 

Q  <  -1.0 

compressed  solid  states 

3.     OVERVIEW  OF  THE  PROGRAMS 
The  decks  available  and  described  in  this  report  are  listed  in  table  2. 

Table  2.     List  of  Available  Program  Decks 


Input 

Output 

p-hydrogen 

e -hydrogen 

P-T 

enthalpy- 

PTENTH 

ETENTH 

entropy 

TPENTR 

ETENTR 

density 

PTDENS 

ETDENS 

Cp 

PTCP 

ETCP 

Cv 

entry  PTCV 

entry  ETCV 

gamma 

entry  PTGAMM 

entry  ETGAMM 

sound  velocity 

PTSOUN 

ETSOUN 

viscosity 

PTVISC 

ETVISC 

thermal  conductivity- 

PTCOND 

ETCOND 

heat  transfer  parameter 

PTLFAC 

ETLFAC 

P-H 

temperature,  and  Q 

PHTEMP 

EHTEMP 

entropy 

PHENTR 

EHENTR 

density 

PHDENS 

EHDENS 

CP 

PHCP 

EHCP 

cv 

entry  PHCV 

entry  EHCV 

gamma 

entry  PHGAMM 

entry  EHGAMM 

sound  velocity 

PHSOUN 

EHSOUN 

viscosity 

PHVISC 

EHVISC 

thermal  conductivity- 

PHCOND 

EHCOND 

heat  transfer  parameter 

PHLFAC 

EHLFAC 

Saturat 

ion  Routines 

P 
vap 

rp                              TT                              TT 

-■-vap'      gas'   "liquid 

PSATH 

EPSATH 

T             1 
vap 

Pvap'   Hgas'   Hliquid 

TSATH 

ETSATH 

Including  the  small  saturation  routines  there  are  36  decks  (the  output  options  C    ,    C    . 

p         v 

and  the  ratio  C    /C     are  combined  into  one  deck).      Each  deck  exists  in  a  CDC  and  an  IBM 
P       v 

version.     Program  listings  are  given  in  Appendix  B. 


3.  1     Method  Used 

The  conceptual  idea  of  the  computer  programs  follows  from  the  use  of  a  T-S  or 
Mollier  diagram.     Using   such  a  chart  to  obtain  property  data  one  visually  or  graphically 
interpolates  between  constant  property  lines.     The  computer  programs   imitate  this  procedure 
on  the  computer.     However,    in  contrast  to  the  visual  procedure  where  any  combination  of 
existing  variables  can  be  used  as  entries  to  the  chart,    the  computer  programs  are  restricted 
to  two  entry  possibilities,  namely  pressure- temperature  and  pressure- enthalpy. 

Each  property  forms  a  surface  with  the  two  entry  variables.     A  pre-determined 
number  of  points  on  this  surface  are  calculated  and  form  the  array  of  data.     The  surface  is 
divided  into  regions,    the  boundaries  of  which  depend  upon  the  relative  linearity  of  the  partic- 
ular property  surface.     The  P,    T,    or  H  spacings  vary  from  region  to  region  to  compensate 
for  the  non-linearity  of  the  property  surface.     In  general,    spacings  will  be  small  where  the 
property  is  varying  quite  rapidly  and  they  will  be  large  where  the  property  is  varying  more 
slowly.     Each  program  consists  of  an  array  of  data  and  an  interpolation  package.     The 
interpolation  package  finds  the  four  points  which  bracket  the  entry  values  and  then  performs 
a  two  dimensional  linear  interpolation  to  obtain  the  value  of  the  property  to  be  returned. 

The  arrays  of  data  are  necessarily  limited.     When  coupled  with  linear  interpolation, 
errors  in  the  values  returned  are  possible.     We  have,    therefore,   determined  the  error 
envelope  for  each  program  and  these  are  displayed  as  departure  graphs. 

3.2  Peculiarities  of  the  P- T  Programs 

One  prominent  feature  shown  in  the  appropriate  phase  diagram,   figure  3,    is  the  vapor 
pressure  curve.      In  all  of  the  P-T  programs  this  line   is  mathematically  thin,    that  is,    the 
property  value  returned  corresponds  to  either  the  liquid  or  the  vapor   side.      In  following  an 
isobar,   an  isotherm,   or  any  other  path  that  crosses  the  vapor  pressure  curve,    one  will 
encounter  sizeable  changes  in  the  different  variables.     In  addition,  for  many  variables  these 
large  changes  also  exist  for  conditions  near  or  just  above  the  critical  point.     Re- stating  the 
peculiarity  a  different  way:     it  is  impossible  to  pick  a  P-T  input  which  will  place  the  coordi- 
nates inside  a  two  phase  region.     Since  it  is  often  necessary  to  obtain  property  values  for 
conditions  of  saturation,    saturated  liquid,    saturated  vapor,   or  both,   we  have  provided  a  way 
of  obtaining  these  values.     The  saturation  routines  listed  in  table  2  are  designed  to  fill  this 
need.     They  are  described  in  Section  3.4. 

3.3  Peculiarities  of  the  P-H  Programs 

In  contrast  to  the  P-T  programs  it  is  quite  possible  to  choose  input  values  of  P  and  H 
which  place  a  point  in  question  into  a  two  phase  region.     From  figure  4  it  is  seen  that  the 
boundaries  of  the  programs  are  set  specifically  to  omit  returns  from  the  liquid- solid  (slush) 
two  phase  region.     Returns  from  the  liquid-vapor  two  phase  region  are,   however,   possible. 
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Three  distinct  possibilities  or  problems  connected  with  these  returns  are  discussed  in  some 
detail  below. 

The  Phase  Designator,    Q,    from  PHTEMP  -  One  problem  often  encountered  is  that  we  need 
to  know  where  we  are  on  the  phase  diagram.      To  solve  this  problem  the  program  PHTEMP 
(EHTEMP)  is  structured  slightly  different  than  all  of  the  other  programs.      The  other  P-H 
programs  carry  the  two  entry  variables  P  and  H  in  the  call  list.      PHTEMP  carries  an 
additional  variable,    O,    in  its  call  list.     A  call  to  PHTEMP  (P,  H,  Q)  and  a  check  of  the  O 
that  is  returned  against  table  1  will  locate  the  point  on  the  phase  diagram.      For  example, 
if  0  <  Q  <  1,    then  the  P-H  input  is  in  the  liquid- vapor  coexistence  region.      For  this  region 
Q  corresponds  to  the  conventional  quality  and  is  defined  as  the  percent  (by  mass)  of  the 
vapor  phase  present,    or 

PQ    =     (H  -  H.  )  /   (H  -  H.      .,) 

liquid  gas  liquid 

where  H, .       .  ,  and  H  are  the  enthalpy  values  at  the  respective  phase  boundary  appropriate 

liquid  gas 

to  the  entry  pressure  (see  also  fig.    4).     As  noted  before,    values  of  Q  for  other  parts  of  the 
phase  diagram  are  given  in  table  1. 

Returns  from  the  Liquid-Vapor  Two  Phase  Region  -  Assume  that  the  input  values  of  P  and  H 
place  the  point  in  question  into  the  liquid-vapor  two  phase  region.     For  this  region  the  entry 
enthalpy  H  is  used  in  all  PH  (EH)  programs  to  find  t,he  phase  designator  or  quality  as 
defined  above,  and  the  phase  designator  or  quality  is  then  used  to  define  the  return,   RE. 
For  temperature,   density,  and  entropy,   these  returns  have  the  conventional  meaning  which 
is  given  in  table  3.     Other  variables,   that  is  k,    T),  tu,   C    ,   and  L,  are  not  defined  for  this 
region  in  the  sense  of  thermodynamics.     However,   the  engineer  often  uses  these  properties 
in  this  region,   and  then  defines  them  according  to  the  specific  application.     Therefore,    the 
returns  given  are  artificial  involving  Q  and  property  values  corresponding  to  the  saturated 
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liquid  and  saturated  vapor  as  defined  in  table  3.     C       and    y    are  infinite  in  this  two  phase 
region;  thus  the  largest  value  available  to  the  machine  is  returned.     The  returns  are 
summarized  in  table  3. 

Property  Values  on  the  Saturation  Boundaries  -  Property  values  of  saturation  are  required  in 
many  applications.     From  the  preceding  discussions  it  is  clear  that  these  values  are  simply 
a  special  case  of  returns  from  the  liquid-vapor  two  phase  region.     A  property  value  on  a 
phase  boundary  is  obtained  from  the  appropriate  PH  (EH)  program  by  using  pressure  and  the 
saturation  enthalpy,   liquid  or  gas,   as  entries.     The  saturation  enthalpies  can  be  obtained 
from  the  two  small  saturation  routines. 

3.4    The  Saturation  Routines 

For  points  on  the  saturation  boundary,    i.  e.  ,    on  the  vapor  pressure  curve  of  figure  3, 
either  P  or  T  is  sufficient  to  specify  the  conditions  of  the  point  in  question  (phase  rule). 
Obviously,     the  pressures  have  to  be  between  1.022  and  187.506  psia  and  the  temperatures 
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Table  3.     Returns  from  the  Liquid-Vapor  Two  Phase  Region 


Prog 

ram 

return  mode 

description  of  return 

PHTEMP 

EHTEMP 

conventional 

RE  =  the  temperature  corresponding 
to  the  input  pressure  (vapor  pressure). 
Additional  output:    the  phase  designator, 
D. 

PHDENS 

EHDENS 

conventional 

RE  =  average  density  of  mixture 

=   1.0/[Q/p               +  (l-Q)/p..      .,] 
L          vapor                        liquid-1 

PHENTR 

EHENTR 

conventional 

RE  proportional  to  quantities  of 

PHCOND 

EHCOND 

artificial 

liquid  and  vapor  present  according 

PHVISC 

EHVISC 

artificial 

to  the  equation 

PHSOUN 

EHSOUN 

artificial 

PHCV 
PHLFAC 

EHCV 
EHLFAC 

artificial 
artificial 

RE  =  Q(tt               -  tt             )  +  tt 

vapor          liquid           liquid 

where    tt     stands  for  the  property 
in  question 

PHCP 

EHCP 

infinite 

RE  =  the  largest  value  available  to 

PHGAMM 

EHGAMM 

infinite 

the  machine 

between  24.845  and  59.356°R.     The  saturation  routines  need  to  be  entered  with  only  one 

variable;  the  corresponding  other  variable  on  the  vapor  pressure  curve  is  returned.     In 

addition,   the  enthalpies  of  the  saturated  liquid,  H, .       .  , ,  and  the  saturated  vapor,   H         ,   are 

liquid  gas 

returned  corresponding  to  the  two  phase  boundaries  of  figure  4.     Saturation  values  for  all 
other  properties  can  then  be  obtained  by  using  pressure  and  the  appropriate  saturation 
enthalpy  as  input  to  the  desired  P-H  program,   as  shown  in  the  example  on  page  16. 


4.      PROGRAMMING  INFORMATION 
4.  1     Deck  Description 

All  of  the  decks  are  coded  as  FORTRAN  functions.      The  decks  are  completely 
self-contained,    including  the  data  arrays  which  are  used  for  interpolation.     Each  deck  may 
be  used  separately  or  in  conjunction  with  others. 

A  prefix  PT  (ET)  on  a  subprogram  name  means  that  it  requires  pressure-temperature 
as  input  variables  while  the  prefix  PH  (EH)  indicates  that  the  input  variables  are  pressure- 
enthalpy.      The  remainder  of  the  deck  name  designates  the  variable  returned.      The  names  of 
the  decks  and  units  for  input  and  output  are  given  in  table  4  for  p-hydrogen  and  in  table  5  for 
e-hydrogen.     Also  shown  in  these  tables  are  data  loading  routines,    storage  requirements, 
and  number  of  cards  in  each  deck. 
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Table  4.     Deck  Description,    p-hydrogen 


Input  variables 
and  units 


>ub-program 
name 

Output  variable 
and  unit 

Extra   sub- 
routine for 
data  loading 

Storage 
requirement 

Number 
of  cards 
in  deck 

PTENTH 

enthalpy 
BTU/lb 

CDC   -  2 

1717 

CDC  258 

IBM  256 

PTDENS 

density 
lb  /ft3 

CDC   -  3 

1460 

CDC  234 
IBM  227 

PTENTR 

entropy 
BTU/lb°R 

CDC   -   1 

1274 

CDC  208 
IBM  208 

PTCOND 

thermal 
conductivity 
BTU  ft/h     ft2  °R 

CDC  -   1 

1156 

CDC  184 
IBM  182 

PTVISC 

viscosity 
lb  h    /ft2 

CDC  -   1 

978 

CDC   197 
IBM  195 

PTSOUN 

velocity  of 
sound  ft/s 

CDC  -   1 

980 

CDC   164 

IBM  165 

PTCP 

PTCV 
PTGAMM 

C    ,    BTU/lb°R 
CP,    BTU/lb°R 
Y  *      NONE 

CDC  -  4 
IBM  -   1 

2151 

CDC  292 

IBM  291 

PTLFAC 

ko.s   co.4/no.4 

CDC   -   1 

1157 

CDC  183 

pressure  psia 

and 
temperature 

•R 


IBM  180 


pressure  psia 

and 
enthalpy 
BTU/lb 
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PHTEMP 


PHDENS 


temperature 

°R 
quality  see  text 

density 
lb /ft3 


CDC   -   1 


917 


1142 


CDC  160 

IBM  161 


CDC   193 
IBM  194 


11 


PHENTR 


entropy 
BTU/lb°R 


CDC   -   1 


1130 


CDC   181 
IBM  180 


12 


PHCOND 


thermal 
conductivity 
BTU  ft/h      ft2°R 


CDC   -   1 


1304 


CDC  213 

IBM  209 


13 


PHVISC 


viscosity 
lb  h    /ft2 


CDC  -   1 


1402 


CDC  266 

IBM  263 


14 


PHSOUN 


velocity  of 
sound,    ft/s 


NONE 


795 


CDC  145 

IBM  148 


15 


PHCP 
PHCV 
PHGAMM 


C  ,  BTU/lb°R 
CP,  BTU/lb°R 
Y,V      NONE 


CDC  -  4 

IBM  -   1 


2077 


CDC  294 
IBM  289 


16 


PHLFAC 


<0-6   C°-4/Tl0-4 
P 


CDC   -   1 


1304 


CDC  215 
IBM  210 


temperature 

•R 


17 


TSATH 


vapor  pressure  NONE 

Hgas,    BTU/lb 
HUq,    BTU/lb 


161 


CDC  29 
IBM  32 


pressure 
psia 


18  PSATH  sat.    temperature         NONE 

Hgas>    BTU/lb 
HUq,    BTU/lb 

13 


161 


CDC  29 
IBM  32 


Table  5.     Deck  Description,    e-hydrogen 


Input  variables 
and  units 


Deck  Sub-program 

#  name 


Output  variable  Extra  sub-  Storage  Number 

and  unit  routine  for  requirement  of  cards 

data  loading  in  deck 


pressure  psia 

and 
temperature 

'R 


ETENTH 


ETDENS 


enthalpy 
BTU/lb 


density 
lb /ft3 


CDC   -   1 


CDC   -  2 


1717 


1460 


CDC  279 

IBM  255 

CDC  229 
IBM  225 


ETENTR 


entropy 
BTU/lb°R 


CDC   -   1 


1274 


CDC  209 
IBM  206 


ETCOND 


thermal 
conductivity 
BTU  ft/h     ft2°R 


CDC   -   1 


1156 


CDC  184 
IBM  182 


ETVISC 


viscosity 
lb  h   /ft2 


CDC   -   1 


978 


CDC  197 
IBM  195 


ETSOUN 


velocity  of 
sound  ft/s 


CDC  -   1 


CDC  164 

IBM  162 


ETCP 
ETCV 
ETGAMM 


Cp,  BTU/lb°R 
Cv,  BTU/lb°R 
Y,       NONE 


CDC   -  3 
IBM  -  1 


2151 


CDC  291 
IBM  293 


ETLFAC 


k°-s  c0-4/"0-4 


CDC   -   1 


1157 


CDC  191 
IBM  180 


pressure  psia 

and 
enthalpy 
BTU/lb 


EHTEMP 


temperature 

•R 
quality  see  text 


CDC   -   1 


917 


CDC  161 

IBM  160 
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EHDENS 


density 
lb /ft3 


CDC  -   1 


1142 


CDC   195 
IBM  194 


11 


EHENTR 


entropy 
BTU/lb°R 


CDC  -   1 


1130 


CDC  182 

IBM  180 


12 


EHCOND 


thermal 
conductivity 
BTU  ft/h     ft2oR 


CDC  -   1 


1304 


CDC  213 
IBM  211 


13 


EHVISC 


viscosity 
lb  h    /ft2 


CDC  -   1 


1402 


CDC  267 
IBM  267 


14 


EHSOUN 


velocity  of 
sound,    ft/s 


NONE 


CDC  141 

IBM  141 


15 


EHCP 
EHCV 
EHGAMM 


C  BTU/lb°R 
Cv,  BTU/lb°R 
Y,       NONE 


CDC  -  2 

IBM  -   1 


2077 


CDC  289 
IBM  291 


16 


EHLFAC 


k0.6    C0.4/T10.4 

p 


CDC   -   1 


1304 


CDC  227 

IBM  210 


temperature 
°R 


17 


ETSATH 


vapor  pressure 
Hgas,    BTU/lb 
Hliq,    BTU/lb 


NONE 


161 


CDC  30 
IBM  30 


pressure 
psia 


18 


EPSATH 


sat.    temperature 
Hgas-    BTU/lb 
Hiiq.    BTU/lb 


NONE 


161 


CDC  30 

IBM  30 
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Since  nearly  every  computer  has  available  some  type  of  FORTRAN  processing 
language,    the  programs  were  written  in  FORTRAN  IV.      Program  listings  for  both  the  CDC 
and  the  IBM  versions  of  all  decks  are  given  in  Appendix  B. 

4.  2     Calling  Sequence 
Since  the  decks  are  individual  FORTRAN  functions  the  calling  sequences  are 
predictably  simple. 

P  =  x.  x 
T  =  y.  y 
H  =   z.  z 
and  then 

RHO  =  PTDENS  (P,  T), 
or  RHOE  =  EHDENS  (P,  H). 

Provided  that  the  values  of  P,    T,    and  H  are  valid,    that  is,    within  the  program  range, 
the  density  of  para  hydrogen  is  stored  into  RHO  or  the  density  of  equilibrium  hydrogen 
into  RHOE. 

The  only  variation  on  this  theme  occurs  for  the  routine  PHTEMP  (EHTEMP)  and 
the  saturation  routines.      The  change  occurs  not  in  the  calling  sequence  but  rather  in  what 
information  is  available  after  the  call  is  made.      Examples  are  given  below. 
P  =  x.  x 
H  =  z.  z 
TEMP  =  PHTEMPfP,  H,  Q)  or  TEMP  =  EHTEMP(P,  H,  Q). 

The  temperature  corresponding  to  the  input  pressure  is  found  in  TEMP.     In  addition, 
the  phase  designator  Q  is  available  through  the  parameter  list,    for  example, 

QUAL  =  Q. 
Note  that  we  did  not  set  Q  prior  to  the  call.      The  saturation  routines  use  very  similar  logic. 

P  =  x.  x 
TVAP  =  PSATH(P,  HG,  HL)  or  TVAP  =  EPSATH(P,  HG,  HL) 

The  temperature  corresponding  to  the  input  pressure  (vapor  pressure)  is  found  in  TVAP. 
In  addition,    the  enthalpy  of  the  saturated  liquid  and  the  enthalpy  of  the  saturated  vapor  are 
available  through  the  parameter  list;    thus 
HGAS  =  HG,       and 
HLIQ  =  HL  . 
The  sequence  with  the  other  saturation  routine  is  as  follows: 

T  =  y.  y 
PVAP  =  TSATH(T,  HG,  HL)  or  PVAP  =  ETSATH(T,  HG,  HL). 

Here  the  vapor  pressure  is  returned  in  PVAP,    and  again  the  enthalpies  of  saturated 
liquid  and  vapor  are  available  through  the  parameter  list.     Next,    assume  the  call  above 
has  been  made,    but  what  really  is  required  is  the  entropy  of  the  saturated  vapor.      The 
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sequence  is  completed  by 

SGAS  =  PHENTR  (PVAP,  HG)  or  SGASE  =  EHENTR  (PVAP,  HG). 

4.  3    Exceeding  the  Range  of  the  Programs 

It  is,    of  course,    easy  to  violate  the  range  of  the  programs  when  setting  the  values 
of  the  input  variables.      The  problem  seems  to  occur  frequently  if  unit  conversions  are 
involved.     We  have  omitted  error  messages  for  cases  of  boundary  violations  because  each 
installation  handles  printed  error  messages  differently. 

For  the  PT  (ET)  and  PH  (EH)  programs,    the  out  of  bounds  variable  is  automatically 
set  to  a  value  acceptable  to  the  program,    and  computations  proceed.      The  value  returned 
is,    of  course,    in  error.     Possible  boundary  violations  are  shown  in  the  phase  diagrams, 
figures  3  and  4,    by  dots,    and  the  arrows  indicate  what  conditions  are  actually  processed. 
A  complete  description  is  as  follows: 

Pressure  P  out  of  bounds 
If  the  entry  pressure  is  less  than  1  psia,    the  value  of  pressure  is  set  to  1. 
If  the  entry  pressure  is  greater  than  5000  psia,    the  returned  value  is  extrapolated  from 
within  the  defined  range. 

Temperature  T  out  of  bounds 
If  the  entry  temperature  is  less  than  the  temperature  appropriate  to  melting  conditions,    the 
temperature  is  set  to  the  melting  temperature  at  the  same  pressure. 
If  the  entry  temperature  is  greater  than  6000°R,    the  value  of  temperature  is  set  to  6000° R. 

Enthalpy  H  out  of  bounds 
If  the  entry  enthalpy  is  less  than  the  enthalpy  appropriate  to  melting  conditions,    the  enthalpy 
is  set  to  the  melting  enthalpy. 

If  the  entry  enthalpy  is  greater  than  25,  000  BTU/lb,    the  value  of  the  enthalpy  is  set  to 
25,  000  BTU/lb. 

Saturation  routines 
If  either  program  is  entered  with  a  value  outside  the  pressure  or  temperature  range  of  the 
vapor  pressure  curve,    the  variable  is  set  to  the  nearest  end  point  of  the  curve  (see  figure  3). 
4.  4    Conversion  to  Other  Systems  or  Installations 

The  programs  as  listed  in  Appendix  B  were  written  to  avoid  known  incompatibilities. 
We  have  tested  them  on  a  CDC  3600,    a  CDC  3800,    and  an  IBM  360.     Other  installations  have 
run  them  on  an  1108  Univac.      Thus,    conversion  requirements  should  be  minimal.     In- 
compatibilities that  we  are  aware  of  include 

DATA  / /  vs.  DATA  ( ) 

ABSF  ( )  vs.  ABS  ( ) 

EXPF  ( )  vs.  EXP  ( ) 

Another  source  of  incompatibilities  arises  from  the  use  of  different  key  board 
symbols  on  certain  types  of  card  punches,    i.  e.  ,  (       )      +      =  vs.  %      <      &  ={  J= 
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The  programmer  should  note  that  the  data  arrays  are  entered  via  DATA  statements. 
The  upper  limit  on  the  number  of  such  statements  varies  from  installation  to  installation. 
The  least  number  that  we  are  aware  of  has  been  used.     Since  this  means  that  in  many  cases 
not  all  of  the  data  can  be  loaded  with  the  actual  program  deck,    we  had  recourse  to  data 
loading  subroutines  which  are  listed  in  tables  4  and  5.     If  a  limitation  is  encountered,    a 
possible  solution  is  as  follows  (Goldberg,    1967): 
Problem:  File  length  not  sufficient  to  include  all  DATA  statements  of  routine 

PTENTH. 
Solution:  Two  segments  used,    1.     BLOCK  DATA,     2.      Function  PTENTH  modified. 

Details:  1.     BLOCK  DATA  (new  for  PTENTH  function)  COMMON/BLKPTH/AA.  .  . 

AO,    followed  by  the  DATA  statements  for  arrays  AA  through  AO  (total 

1183  values  equivalent  to  H  array). 

2.     Main  function  PTENTH  is  revised  to  include  all  DATA  arrays 

except  H  (1183).      The  H  array  is  linked  to  the  main  function  through 

named  common:     COMMON/ BLKPTH/H(1 183). 

4.  5     Pitfalls,    Common  Errors  in  Applications 

The  earlier  version  of  these  decks  has  been  available  for  nearly  five  years. 

Several  points,    where  users  have  encountered  difficulties,    are  described  below: 

Storage  Limitation  -  Each  program  is  designed  to  be  used  independently  of  all  others. 

Specify  only  those  decks  actually  required.      Total  storage  required  is  specified  in  tables 

4  and  5.     Storage  requirements  range  from  795  to  2151. 

Program  Range  -  Observe  the  limits  of  the  programs  1  -  5000  psia,    24.  845  -  6000°R, 

-130  -  25,  000  BTU/lb.     Valid  answers  will  not  be  obtained  if  these  bounds  are  violated. 

Liquid-Vapor  Two  Phase  Conditions  -   There  is  no  built-in  warning  in  any  of  the  PH(EH) 

decks  telling  the  user  that  he  is  in  the  liquid-vapor  two  phase  region.      This  means  the  user 

should  check  the  value  of  Q  to  confirm  the  conditions  he  is  trying  to  process. 

Liquid-Vapor  Two  Phase  Returns   -   Returns  from  this  region  for  the  variables  k,    T|,   cju,    C    , 
p 

C    ,  y,  and  L,  are  artificial,  as  defined  in  table  3.     The  programs  will  return  numbers  for 
these  variables,   but  what  usually  is  desired  are  not  these  returns  but  rather  the  property 
values  for  saturated  liquid  or  saturated  vapor. 

Errors  -  The  two  categories  of  errors,    interpolation  error  and  estimated  source  data  errors, 
are  discussed  in  detail  in  Section  6.     The  user  has  to  be  aware  of  the  limitations  inherent 
in  linear  interpolation.     The  detailed  graphs  of  maximum  interpolation  errors  given  in 
Appendix  C  should  alert  the  user  to  possible  problems,   and  serve  as  an  aid  in  solving  those 
problems. 

4.  6    Debugging  Aids 
In  addition  to  the  program  listings  given  in  Appendix  B,   we  have  found  that  a  set  of 
selected  values  for  the  various  properties  is  indispensable.      In  Appendix  A  we  list  values 
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obtained  from  the  PT(ET)  programs.     Values  for  all  properties  for  a  representative  set 
of  temperatures  are  printed  for  isobars  of  1,    15,    50,    150,    500,    1500,    and  4500  psia. 
These  skeleton  tables  are  to  be  used  to  establish  approximate  values,    to  make  an  order  of 
magnitude  check,    to  check  units,    and  to  debug  the  running  of  the  programs  at  other 
installations. 

5.      STRUCTURAL  DETAILS  OF  THE  PROGRAMS 

As  mentioned  before,    each  program  can  be  visualized  as  an  array  of  data  coupled  with 
an  interpolation  scheme.     Each  property  value  entered  into  the  data  array  corresponds  to  a 
particular  pair  of  P-T  or  P-H.      The  spacing  or  increment  for  each  entry  variable  is 
independent  of  the  other  and  varies  over  the  surface  of  the  property.      The  total  network  of 
spacings  is  called  the  grid  and  a  sequence  of  equal  spacings  defines  a  region. 

5.  1     Grid  and  Region  Layout 

The  grid  for  each  variable  consists  of  12  to  25  different  regions.      The  size  of  each 
region  and  its  relation  to  neighboring  regions  is  determined  by  the  desired  interpolation 
error  and  by  the  relative  linearity  of  the  particular  property  surface.     A  plot  of  constant 
property  lines  is  used  as  an  aid  in  layout.     An  example  of  grid-region  layout  for  p-hydrogen 
is  shown  in  figure  5  where  lines  of  constant    C       (in  BTU/lb-°R)  are  plotted  for  entry 
variables  P  and  H.      Lines  of  constant    C       and  regions  G  through  Q  are  shown.     Within  each 
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region  increments  of  the  entry  variables,   P-T  or  P-H,  are  fixed.     The  spacings  used  for 

the  various  regions  were  intended  to  yield  a  1%  interpolation  error;  however,    subsequent 

checks  show  that  we  have  not  always  achieved  this  goal.     A  data  array  developed  according 

to  this  scheme  will  vary  from  variable  to  variable;  it  will  minimize  the  number  of  points 

required  to  represent  the  surface  of  the  variable  to  be  interpolated. 

While  in  general  the  property  values  entered  in  the  data  arrays  are  exact,    there  are 

a  few  cases  where  these  entries  have  been  adjusted  to  yield  better  results  for  the  inter- 
polation.    One  example  is  entries  at  low  pressures,    that  is,    some  regions  are  extended  to 
negative  pressures  and  the  entries  at  these  negative  pressures  are  bogus  values.      They  do, 
however,    insure  a  more  precise  interpolation  down  to  the  lower  limit  of  pressure.     Region 
B  shown  in  the  partial  grid  layout  of  figure  6  is  a  typical  example.      The  lower  pressure  limit 
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for  interpolation  from  this  region  is  80  psia  and  entries  in  the  data  array  at  -100  psia  are 
bogus  values  (shown  as    ©).     A   similar  technique  better  described  as  a  partial  overlap  of 
region  boundaries  is  used  in  Region  A,   C,   and  D.     For  best  accuracy  values  at  15  psia  are 
interpolated  from  Region  D,   not  from  Regions  A,   B,   or  C.     Region  E  is  the  only  one 
illustrated  in  which  bogus  values  are  not  used.     The  entire  sequence  illustrates  how  difficult 
it  is  to  follow  the  large  variation  in  enthalpy  in  the  region  of  dissociation. 

Interpolation  along  the  two  phase  boundaries  also  involves  the  use  of  bogus  points 
that  have  no  counter  part  on  the  surface  of  the  property  in  question.     For  the  P-T  programs, 
the  points  are  across  the  two  phase  boundary  while  for  the  P-H  programs  they    have  been 
added  inside  the  two  phase  boundaries.     In  either  case,   these  bogus  points  are  used  only  as 
interpolation  aids  to  get  from  the  single  phase  region  right  up  to  the  two  phase  boundary. 

5.  2    Speed 

The  programs  gain  their  operational  speed  from  performing  two  dimensional  linear 
interpolation  in  a  fixed  array  of  data.     Actual  speed  of  computation  depends,    of  course,    on 
the  particular  installation.     At  our  installation,    a  CDC  3800,    returns  run  from  30,  000  to 
60,  000  per  minute.      The  data  used  for  interpolation  have  been  placed  into  compact  DATA 
statements  to  provide  rapid  access.     Access  to  any  given  region  is  just  as  fast  as  to  any 
other  region  because  of  the  way  searching  is  arranged.      The  first  branching  divides  the 
number  of  regions  in  half,    the  second  in  half  again,    and  so  on.     Speed  is  reduced  by  about 
a  factor  of  2  if  lookup  of  the  phase  boundaries  is  required  for  the  point  under  consideration. 
The  programs  return  a  single  answer  with  two  variables  (three  for  PHTEMP)  going 
through  the  transfer  vector. 

5.  3    Continuous  Functions 

In  the  approach  taken  here,    continuous  functions  cannot  be  obtained  because  the 
increments  of  the  entry  variables  change  when  crossing  region  boundaries.      The  problem  is 
illustrated  by  two  simple  examples  shown  in  figure  7.     In  case  I,    points  1,  2,  3,  4,    and  5  lie 
on  the  boundary  of  the  lower  region,    while  points  1,  3,    and  5  are  points  on  the  boundary  of 
the  adjacent  upper  region.      Linear  interpolation  along  the  boundary  for  each  region  is  shown 
by  the  straight  lines.     If  points  2  and  4  are  adjusted  to  2'  and  41,    interpolation  in  both  regions 
will  be  the  same  and  jumps  in  the  value  of  the  property  returned  will  not  occur  across  this 
boundary.      In  case  II,    points  1,  2,  4,    and  5  are  in  the  lower  region  while  1,  3,    and  5  lie  in 
the  adjacent  upper  region.      The  values  of  2  and  4  cannot  be  adjusted  to  avoid  a  jump  in  the 
property  interpolated*.     An  actual  example  where  boundaries  cannot  be  made  to  match  can  be 
seen  in  the  transition  from  30  psia,    Region  D,    to  30  psia,    Region  C,    of  figure  6. 


*    We  cannot  adjust  the  value  of  point  3  as  the  interpolation  error  at  3  would  be  larger  than 
the  desired  1%. 
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boundaries  can  be   made  continuous. 

Figure  7.      Matching  of  Values  at  Boundaries 

5.  4    Modifications,    Future  Developments 

Changes  are  a  way  of  life  with  computer  programs.      Improvements  will  be  dictated 

by  the  user  —    they  may  arise  from  better  source  data.      One  possible  modification  is 

described  below. 

Return  of  Values  on  the  Two  Phase  Boundary  -  In  many  applications  involving  the  PH  (EH) 

programs,    some  of  the  returns  are  from  the  two  phase  region.      In  this  case,    the  property 

values  on  the  respective  boundaries  are  required.      Since  the  programs  calculate  the  values 

at  the  boundaries  for  the  regular  return,    a  slight  change  in  program  structure  will  yield 

the  desired  answers.      The  changes  are  illustrated  by  the  following  example  for  entropy: 

FUNCTION  PHENTR  (PRES,  ENTH)  change  to 
FUNCTION  PHENTR  (PRESS,  ENTH,  SGG,  SLL) 
SGG  =  SGAS 
SLL  =  SLIQ 


and  add  at  25  +  2 


In  other  words,    the  new  logic  is  similar  to  that  used  in  PHTEMPjP,  H,  Q)  and  is  accomplished 
at  the  expense  of  having  two  additional  variables  pass   through  the  transfer  vector. 
Future  Possibilities   -  It  is  interesting  to  speculate  on  improvements  that  might  be  accomplish- 
ed in  the  future  based  on  what  we  have  learned  about  the  linear  interpolation  scheme.     An  ad- 
vanced scheme  of  interpolation  would  save  storage,    improve  accuracy,    and  allow  a  large 


combination  of  variables.      It  would  be  based  on  the  following  elements: 

21 


1.  The  development  of  a   single  program  to  include 

a)  an  interpolation  package 

b)  a  data  array  of  fixed  size,    about  350  points  per  property. 

2.  The  use  of  higher  order  or,    perhaps  spline  interpolation  in  the  interpolation 
package. 

3.  The  use  of  the  distribution  in  grid  entries  obtained  during  the  present  error 
analysis  (see  for  example  tables  8  and  9)  to  create  an  optimum  point  distribution 
for  the  new  array  of  data. 

4.  A  separate  program  which  would  call  the  generating  functions  (see  Section  6.  8) 
and  create  the  array  of  data  desired. 

6.     ESTIMATION  OF  ERRORS 

Terms  such  as  accuracy,   precision,  uncertainty,   and  error  are  used  in  scientific  and 
engineering  papers  with  many  shades  of  meaning.     To  avoid  ambiguity,    the  terms  used  in 
this  report  are  defined  as  follows: 

First,   we  define  error  as  the  quantitative  difference  which  exists  between  the  true 
value  and  the  value  which  is  measured,    calculated,   or  otherwise  determined. 

Second,  uncertainty  is  defined  as  an  estimate  of  the  error  that  may  exist  between 
reported  and  true  values.     Since  usually  we  do  not  know  the  true  value,   the  estimate  is  by 
nature  a  variable  quantity.     The  estimate  of  error  or  uncertainty  will  vary  according  to  the 
information  available  and  the  manner  in  which  the  estimate  is  obtained. 

Statistical  methods  normally  are  used  to  arrive  at  estimates  of  error  or  uncertainties. 
However,    statistical  methods  are  valid  only  for  random  processes,   and  therefore  uncertain- 
ties based  on  statistical  methods  are  estimates  based  on  the  presence  of  random  errors  only. 
In  contrast  to  statistics,   the  nature  of  thermophysical  data  is  such  that  systematic  errors 
usually  are  present,    i.e.  ,    errors  which  are  not  random.     Not  only  are  systematic  errors 
present  but  often  they  are  of  much  greater  importance  than  the  truly  random  errors. 

To  be  realistic  in  assigning  uncertainties  we  have  combined  the  estimates  of  both 
random  and  systematic  errors,   except  where  the  uncertainties  are  labeled  "statistical 
uncertainties"  in  which  case  only  random  errors  are  given.     To  minimize  the  risk  of  a  true 
value  being  outside  of  the  interval  specified  by  our  assigned  uncertainty,    the  assignment  of 
uncertainties  is  large  enough  to  give  a  99%  chance  of  including  the  correct  value,    i.e.  ,    a 
3  a  level.     The    a     is  used  here  in  the  conventional  sense  of  standard  deviation. 

Differentiation  is  made  between  uncertainties  in  the  data  of  the  sources  chosen  and 
interpolation  error,    i.e.  ,   the  difference  between  the  values  returned  by  linear  interpolation 
and  the  source  chosen.     Maximum  interpolation  errors  and  estimated  source  errors  for  each 
program  are  displayed  graphically  in  Appendix  C  in  three  ranges  of  temperature:     25  to 
180°R,    180  to  1200°R,   and  1200  to  6000°R.     The  ranges  of  enthalpy  for  the  PH  and  EH  type 
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programs  correspond  closely  to  the  temperature  ranges  just  mentioned.     All  plots  are  given 
in  percentages  except  for  the  low  temperature  range  of  enthalpy,    i.e.  ,    25  to  180°R  for 
program  PTENTH  (ETENTH).     The  enthalpy  in  this  range  of  temperature  can  assume  both 
negative  and  positive  values   (see  also  fig.    4);   therefore,    the  interpolation  and  source  errors 

for  this  range  only  are  plotted  in  BTU/lb.     The  symbol in  the  plots   signifies  that  the 

PVT  surface  for  these  temperatures  is  extrapolated. 

Estimated  source  errors  are  discussed  in  some  detail  in  Section  7.     The  estimated 
source  errors  are  established  by  the  authors  of  the  source  data,   by  comparison  with  experi- 
mental values,    by  noting  the  spread  of  values   if  several  experimental  sets  exist,    or  by 
intercomparison  of  several  methods  of  calculation  if  no  experimental  values  are  available. 
For  the  high  temperature  extrapolations  we  combine  the  estimates  of  Woolley  (1972)  for  the 
ideal  gas,    our  estimates  of  errors  for  the  real  gas  contribution,   and  an  arbitrary  estimate 
for  errors  due  to  dissociation. 

Interpolation  errors  are  apt  to  intrude  on  most  applications.     These  errors  are 

nominally  3%  except  for  the  critical  region  and  the  region  of  dissociation  where  departures 

from  best  values  become  large.     The  critical  parameters  are     P     =   187.  506  psia, 

c 

T     =  59.357°R,   and    H     =   16.  455  BTU/lb  for  p-hydrogen  or    H     =   31 .  225  BTU/lb  for 
c  c  c 

e-hydrogen.     The  onset  of  dissociation  varies  from  2700°R  at  the  lower  pressures  to  5400°R 
at  the  highest  pressures.     The  interpolation  errors  plotted  are  the  envelope  of  maximum 
error  in  interpolation  at  the  temperature  or  enthalpy  in  question.     They  are  established  from 
an  intercomparison  of  values  returned  by  the  interpolation  routine  and  the  generating 
functions.     The  comparisons  are  made  at  a  large  number  of  temperatures  and  a  large  number 
of  pressures  for   each  temperature,    selecting  the  largest  deviation  in  a  given  range  of  tem- 
perature or  enthalpy  for  the  plot.     The  distribution  of  temperatures  and  pressures  at  which 
the  intercompar isons  were  made  are  established  from  the  grids  of  the  properties  by 
collecting  values  at  1/2  of  each  normal  grid  increment  and  eliminating  duplicates.     As  an 
example,    each  program  was  checked  in  the  low  temperature  range  at  every  combination  of 
P  and  T  given  in  tables  6  and  7. 

Table  6.     Low  Temperature  Error  Grid,   Pressures,   psia 

1.022  36.740  146.960  225.000  308.616  418.836  734.800        1763.520 

1.500  44.083  154.308  230.000  315.964  426.184  771.540       2000.000 

2.000  51.436  161.656  235.000  323.312  433.532  808.280       2500.000 

3.000  58.784  169.004  240.000  330.660  440.880  845.020       2645.280 

4.000  66.132  176.352  245.000  338.008  459.250  881.760       3233.120 

5.00J  73.480  180.026  250.000  345.356  477.620  918.500       3820.960 

6.000  80.628  183.700  255.000  352.704  495.990  955.240        4408.800 

7.000  88.176  187.370  260.000  360.052  514.360  991.980       5000.000 

8.000  95.524  190.000  265.000  367.400  532.730  1028.720 

9.000  102.872  195.000  270.000  374.748  551.100  1065.460 

10.000  110.220  200.000  275.000  382.096  569.470  1102.200 

11.022  117.568  205.000  280.000  389-444  587.840  1138.940 

14.696  124.916  210.000  285.000  396.792  624.580  1175-680 

22.044  132.264  215.000  293.920  404.1U0  661.320  1322.640 

29.392  139.612  220.000  301.268  411.488  698.060  1<*69.60G 
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Table  7.     Low  Temperature  Error  Grid,    Temperatures,    °R 

25.J0Q  kO. 700  55.100  62.500  70.500  90.000 

26.000  1+1.600  56.000  63.000  70.571  93.600 

27.500  «+2.50Q  56.500  63.500  71.500  97.200 

29.uJ0  <+3.<+00  57.000  6<+.000  72. GOO  100.800 

30.500  <+<+.300  57.500  6<+.500  72.900  10<+.<+00 

30.830  1+5.200  58.000  65.000  73.800  108.000 

31.700  <+6.1Q0  58.500  65.500  7*+. 700  112.500 

32.6UJ  1+7.000  59.000  66.000  75.600  117.000 

33.500  *+7.900  59.<»00  66.500  76.500  121.500 

3<+.i+Ga  <+8.800  59*625  67.000  77.<*00  126.000 

35.300  <+9.7Q0  59.850  67.500  78.300  135.000 

36.230  50.600  60.000  68.000  79.200  !<+<+. 000 

37.100  51.500  60.500  68.500  80.000  153.000 

38.000  52.<+00  61.000  69.000  82.500  162.000 

38.900  53.300  61.500  69.500  85.000  171.000 

39.800  5<+.2QQ  62.000  70.000  86.<+00  180.000 

For  the  most  part,  users  will  be  concerned  only  with  the  interpolation  error  because 
in  most  cases  the  uncertainties  in  the  source  data  are  negligible  in  comparison  to  the  inter- 
polation error.     Interpolation  errors  can  be  established  quickly  by  visual  estimate  from  the 
appropriate  plot.     If  the  generating  functions  (Section  7.8)  are  used,    only  the  estimated 
source  errors  need  to  be  considered. 


7.      SOURCES  OF  PROPERTIES  DATA 

To  calculate  values  of  the  various  properties,    we  require  the  thermodynamic  functions 
of  the  ideal  gas,    a  description  of  the  PVT  surface,    a  model  to  calculate  the  effects  of 
dissociation,    and  descriptions  of  the  behavior  of  viscosity  and  thermal  conductivity.     The 
major  references  used  are  shown  in  figure  8  and  are  discussed  in  more  detail  below. 

7.  1     Ideal  Gas  Properties 

As  sources  for  the  thermodynamic  properties  of  the  ideal  gas  we  have  considered 

Woolley,    et  al.    (1948), NBS  Circular  564  (Hilsenrath,    et  al.,  1955)  and  NBS  Monograph  20 

(Haar,    et  al.,  1961).     The  first  two  sources  are  virtually  identical;    they  differ  only  because 

Woolley  includes  the  effect  of  nuclear  spin  in  his  tabulations  for  entropy  (i.  e.  ,    2  R  In  2). 

Normally  the  nuclear  spin  contributions  are  omitted;  both  Circular  564  and  Monograph  20  are 

tabulated  according  to  this  convention.     Slight  differences  in  thermal  properties  between 

Circular  564  and  Monograph  20  exist.     They  can  be  detected  at  temperatures  above  1260°R 

(700  K);  they  increase  as  the  temperature  increases.     For  C°     the  difference  is  about  1%  at 

P 
5000°R  (2800  K)  while  for  H°  -  Eg     the  difference  is  44.6  BTU/lb  (209  J/mol).     A  more 

recent  calculation  (Baehr,    et  al.,  1966)  is  very  close  to  the  earlier  values  of  Woolley, 
et  al.    (1948).     We  have,    therefore,    elected  to  use  the  properties  given  by  Woolley  for  para- 
hydrogen  up  to  3600°R   (2000  K).     At  this  temperature  Woolley  discontinues  the  entries  for 
parahydrogen  but  continues  the  entries  for  normal  hydrogen  on  up  to  9000°R  (5000  K).     At 
3600°R,   differences  between  the  ideal  gas  properties  of  para  and  normal  hydrogen  are 
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SOURCE  DATA  -  IDEAL   GAS  PROPERTIES 
Woolley,  et.  a/.,  NBS  J.  Res.  (1948) 
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Figure  8.        Primary  Sources  of  Data 

negligible  except  in  the  case  of  entropy  where  the  difference  is  1.367  BTU/lb-°R.     The 
difference  represents  the  change  in  order  or  configuration  between  para  and  normal  hydrogen. 
Above  3600°R  the  next  temperature  entry  in  Woolley's  tables  is  5400°R  (3000  K).     For  tem- 
peratures of  5400°R  and  above,   the  ideal  gas  properties  of  normal  hydrogen  were  used  for 
both  para  and  equilibrium  hydrogen  since  at  these  temperatures,    equilibrium  and  normal 
hydrogen  properties  are  identical. 

While  the  entries  for  equilibrium  hydrogen  are  perfectly  proper,  for  parahydrogen 
this  means  that  we  assume  complete  conversion  to  normal  hydrogen  between  3600°  and 
5400° R.     A  higher  order  interpolation  insures  that  the  transition  in  entropy  for  the  para- 
hydrogen decks  is  smooth  and  continuous.     It  becomes  necessary  to  switch  to  the  normal 
mixture  at  some  point  because  the  dissociation  calculations  (Section  7.5)  are  based  on 
normal  hydrogen.     It  should  be  noted  that  the  present  programs  are  oriented  toward  a  partic- 
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ular  process  which  begins  in  the  cryogenic  range  and  ends  at  high  temperatures.     It  would 
not  be  proper  to  use  the  parahydrogen  programs  on  a  process  which  begins  in  the  cryogenic 
range,   goes  to  high  temperatures,  and  then  is  returned  to  the  cryogenic  range,  because 
most  likely  the  hydrogen  would  not  be  returned  to  the  para  form  on  cooling. 

7.2    The  PVT  Surface  at  Low  Temperatures 

Values  used  in  the  data  arrays  for  temperatures  up  to  180°R  (100  K),  for  all  prop- 
erties except  viscosity  and  thermal  conductivity,   were  obtained  from  NBS  Monograph  94 
(Roder,    et  al.,  1965).     The  uncertainties  in  source  data  were  taken  from  this  reference  and 
are  displayed  as  estimated  source  errors  in  Appendix  C  for  temperatures  from  25  to  180°R 
(-130  to  425  BTU/lb  or  675  BTU/lb). 

7.3    The  PVT  Surface  at  Intermediate  Temperatures 

Initially,   we  considered  using  the  exponential  equation  of  state  (Woolley,   et  al.,  1948). 
The  difficulties  with  this  equation  are  that  it  cannot  be  integrated  in  closed  form  to  find  the 
derived  properties,    it  has  a  limited  range,    500  to  1200°R  (273  to  673  K),   and  the  uncertain- 
ties in  the  derived  properties  cannot  be  calculated  easily. 

As  an  alternative  we  considered  a  modified  Benedict- Webb-Rubin  (1940)  equation  of 
state.     The  advantages  of  an  equation  of  state  of  this  type  are  that  the  functions  are  internally 
consistent  and  continuous.     In  addition,    statistical  uncertainties  for  both  PVT  and  derived 
properties  can  be  obtained  from  the  least  squares  estimation  of  the  parameters.     Various 
modifications  of  this  equation  have  been  used  to  correlate  the  PVT  surfaces  of  nitrogen 
(Strobridge,    1962),   helium  (Mann,    1962),    carbon  monoxide  (Hust  and  Stewart,    1963),    neon 
(McCarty  and  Stewart,    1965),   argon  (Gosman,   et  al.,  1969),  and  oxygen  (Stewart,    1966).     The 
equation  has  also  been  used  earlier  on  parahydrogen  (Roder  and  Goodwin,    1961)  and  on  an 
intercomparison  of  normal  and  parahydrogen  (Hust  and  Stewart,    1965). 

The  PVT  data  used  to  estimate  the  parameters  of  the  equation  of  state  included  the 
parahydrogen  isotherms  from  50  to  100  K  (Goodwin,   et  al.,  1963),   all  of  the  PVT  sources 
surveyed  by  Woolley,    et  al.    (1948),   and  the  isotherms  of  normal  hydrogen  reported  by 
Michels,   et  al.    (1959).     The  equation  is  as  follows: 


P    =     pRT  +  ps    (N,T  +  Ns+  N3/T  +  N4/T3  +  N5/T4) 

+  p3    (N6T3  +  N7T  +  N8)  +  p3   (N9/T2  +  N,  0  /T3  +  Nxl/T4)  e    1?P 
+  N18    p4  T+  p5    (N13/T3  +N14/T3  +N15/T4)  eN^  P 
+  N,  e    p8    . 

The  coefficients  are  given  in  table  6  for  P  in  atm,   p   in  mol/4,   T  in  K,   and  R  =  0.0820535 

^-atm/mol-K. 
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Table  8.     Coefficients  of  the  Equation  of  State 


Nx  =  1.  1389049685  x    1CT3 

Ns  =       1.8072093722  x  10"1 

N3  =  -5.  309.7164419  x  101 

N4  =  2.3690885344  x  103 

N5  =  -1.6306395805  x  10s 

Ne  =       3.2453595439  x  10_s 

N7  =  -6.3275640592  x  10~7 

N8  =       2.3115255859  x  10~3 

NQ  =       3.4448764044  x  101 


N. 


l  o 


N, 


N 


l  a 


1  3 


nl  4 

Nib 
N, 


Nis 


N 


1  7 


=     -3. 7301781349  x  103 
1.0789473341  x  105 
8.0475741674  x  10-7 
1.6581404268  x  10_1 
=    -1. 8158230417  x  101 
=       4. 9623738040  x  102 
=       5.262288563    x    10-8 
=       -.0018 


The  equation  is  valid  over  the  range  of  the  data  and  was  used  from  180  to  1200°R. 
Unfortunately  the  major  sources  disagree  considerably  at  the  common  boundary  of  180°  R 
(see,  for  example,  Roder,   et  al.,  1965)  and  the  resulting  jumps  in  the  properties  are  seen 
clearly  in  the  plots  of  maximum  interpolation  error  of  Appendix  C.     A  comparison  of 
calculated  statistical  uncertainties  for  PVT  and  derived  properties  with  the  results  of  Roder, 
et  al.    (1965),   Woolley,   et  al.    (1948),   and  Michels,    et  al.    (1959,    1963)  is  shown  in  figure  9 
where  the  base  line  is  the  present  equation  of  state.     The  comparison  confirms  that  the  values 
obtained  for    o    are  realistic  for  both  PVT  and  derived  properties.     For  temperatures  from 
180  to  1200°R  and  enthalpies  from  425  to  5000  BTU/lb  the  calculated  uncertainties  are  shown 
as  source  errors  in  the  plots  of  Appendix  C. 

7.4.     The  PVT  Extrapolation 

The  highest  temperature  for  experimental  PVT  values   is  673  K.     This  means  the  PVT 
surface  from  1200  to  6000°R  has  to  be  extrapolated.     At  these  temperatures  the  gas   is 
nearly  ideal,   that  is,    only  the  lower  virial  coefficients  contribute  to  the  equation  of  state. 
As  a  first  step  one  has  to  know  how  the  virial  coefficient  B  varies  as  a  function  of  temper- 
ature,  and  then  how  much  the  B(T)  contributes  to  the  equation  of  state.     The  B  for  hydrogen 
attains  its  maximum  value  near  300  K.     Based  on  what  is  known  for  helium  (Keesom,    1942), 
the  B  should  decrease  at  higher  temperatures.     Woolley,    et  al.    (1948)  based  the  temperature 
dependence  of  B  on  this  behavior  and  predicted  the  actual  temperature  dependence  from  the 
Lennard- Jones   (12-6)   intermolecular  potential  function.     Considering  an  extreme  value, 

that  is,    5000  psia   (340  atm)  and  5000°R   (2800  K),    the  exponential  equation  of  state  (Woolley, 

3 
et  al.,  1948)  yields  a  density  of  0.337  lb/ft      (60  Amagat).     In  the  range  in  which  this  equation 

is  valid,    360  to  1200°R   (200  to  673  K),    the  virial  B  contributes  less  than  0.  2%  to  the  density 

for  densities  less  than  60  Amagat.     The  contribution  at  higher  temperatures  can  only  be  less; 

therefore,   an  extrapolation  based  on  an  equation  of  state  which  includes  only  the  second 

virial  coefficient  is  acceptable  in  terms  of  error. 
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Figure  9.     Comparison  of  PVT  and  Derived  Properties 
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Next  we  estimate  the  apparent  error   in  the  extrapolation  of  B(T).     For  5400°R 

3 
(3000  K),    extrapolation  of  the  Woolley,    et  al.    (1948)  calculation  yields   14.55  cm    /mol  for  B 

3 
while  Fisher   (1965)  obtained  12.93  cm    /mol.     Fisher's  calculation  is  based  on  a  "softer" 

intermolecular  potential  function  (i.e.  ,   a  9-6)  and  he  determines  the  molecular  parameters 

from  the  high  temperature  viscosity  measurements  of  Guevara  and  Wageman  (1965).      If 

hydrogen  at  this  temperature  behaved  as  an  ideal  gas,    B  would  be  zero;  thus  a  reasonable 

3 
estimate  for  the  error  in  B  at  5400°R  (3000  K0  is  ±  3  cm    /mol. 

For  the  data  arrays  we  extrapolate  the  equation  of  state   (Section  7.3).      Since  the 

equation  was  fitted  over  the  entire  range  of  data  up  to  1200°R  (6  73  K),    only  one  additional 

constraint  (Hust  and  McCarty,    1967)  was  required  to  yield  acceptable  results.     The  equation 

3  , 
was  constrained  to  yield  14.55  cm    /mol  for  B(T)  at  5400°R  (3000  K).      Z  is  approximately 

1 .  04  at  5400°R  and  5000  psia,    thus  the  real  gas  contribution  is  0.04.     Woolley  (1972) 

estimates  that  the  real  gas  contribution  is  uncertain  by  no  more  than  10%;  therefore,    we 

estimate  the  uncertainty  in  density  to  be  0.8%  for  these  conditions.     Uncertainties   in  the 

derived  properties  are  estimated  in  a  similar  fashion  assuming  an  uncertainty  of  no  more 

than  10%  for  the  real  gas  contribution,    i.e.  ,    the  contribution  from  non  ideality.     This 

estimate  assumes  virtually  no  error   in  the  ideal  gas  thermodynamic  functions  and  virtually 

no  error   in  the  dissociation  calculation  (Woolley,    1972). 

7.5     Dissociation 

At  the  highest  temperatures  considered  in  this  report,    molecular  hydrogen  will 
dissociate  into  atoms.     Dissociation  is  both  temperature  and  pressure  dependent.      To  include 
the  effects  of  dissociation  several  assumptions  are  required. 

1.  Both  ortho-para  and  dissociation  equilibria  occur  instantaneously. 

2.  Dissociation  occurs  according  to  the  equilibrium  constant  for  normal  hydrogen. 

3.  The  property  values  can  be  obtained  by  adding  the  dissociation  effects  calculated  for  the 
ideal  gas  to  the  extrapolated  properties  of  the  real  gas. 

Fortunately,    the  rate  of  ortho-para  conversion  is  quite  high  (Farkas,    1935).     It  is 
proportional  to  the  amount  of  atomic  hydrogen.     The  dissociation  rates  are  apparently  also 
quite  high  (Sutton,    1962).     Assumption  1   is  approximated,   and  thus  assumption  2   is  plausible. 
Hence,    it  is  not  necessary  to  consider  the  "hypothetical"  dissociation  constant  for  para- 
hydrogen  (Farkas,    1935)   even  though  the  hydrogen  is  initially  para. 

The  dissociation  constant  is  calculated  from  the  ideal  gas  free  energies.     As   indicated 
previously,    the  inclusion  or  exclusion  of  nuclear  spin  will  affect  entropy  as  well  as  free 
energy.     The  tables  of  atomic  hydrogen  from  NBS  Circular  564  (Hilsenrath,    et  al.  ,    1955) 
and  the  tables  for  molecular  normal  hydrogen  from  NBS  Monograph  20   (Haar  ,    et  al.  ,    1961) 
were  used.     These  tables  according  to  convention. do  not  include  the  effect  of  nuclear   spin. 
Equation  3.5  for  free  energy  given  by  Woolley,    et  al.,  (1948)   is   equivalent  to  the  table  for 
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atomic  hydrogen  in  Circular  564  but  includes  nuclear  spin.     The  equations  used  to  calculate 
dissociation  effects  of  the  ideal  gas  are  those  given  by  Woolley  (1955). 

The  effects  of  dissociation  on  a  variable  such  as  enthalpy  are  large.     They  are 
illustrated  for  a  limited  range  of  temperatures  and  a  pressure  of  50  psia  in    figure  10  (see 
also  figures  11  and  12).     Enthalpy  differences  are  plotted  vertically  and  the  dissociation 
contribution  ranges  from  2%  to  10%  of  the  total  enthalpy.     Also  shown  is  the  linear  inter- 
polation provided     by  program  PTENTH(P,T). 
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Figure   10.      The  Effect  of  Dissociation  on  Enthalpy 


7.6    Viscosity 
The  primary  experimental  sources  to  consider  are    Diller   (1965)  for  low  temperatures, 
and  Trautz  and  Zink  (1930)  and  Guevara  and  Wageman  (1965)  for  high  temperatures. 
Differences  in  the  latter  two  sources  have  now  been  resolved   in  favor  of  the  more  recent 
data  (Guevara,    et  al.  ,    1971)  allowing  a  critical  evaluation  of  viscosity  and  thermal 
conductivity  coefficients  for  the  dilute  gas   (Hanley,    et  al.  ,    1970).      The  coefficients  are 
represented  by  standard  kinetic  theory  expressions  and  are  based  on  a  model  intermolecular 
potential  function. 
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For  temperatures  up  to  180  °R  for  both  dilute  and  dense  gas,  the  analytical  represen- 
tation of  the  results  by  Diller  (1965)  was  chosen  as  source  data.     At  intermediate  and  high 
temperatures   including  dissociation,    the  evaluation  of    Hanley,    et  al.  ,  (1970)  was  chosen  for 
the  dilute  gas  while  the  dense  gas  contributions  are  calculated  from  the  theory  of  Enskog 
(Hanley,    et  al.  ,    1972).     The  dense  gas  contribution  vanishes  at  1000  K,    that  is,    above 
1000  K  the  viscosity  is   independent  of  pressure  or  density  but  depends  only  on  temperature. 
At  yet  higher  temperatures  with  the  onset  of  dissociation  the  viscosity  becomes  again 
pressure  dependent  as  shown  in    figure   11   taken  from  Hanley,    et  al.    (1970).      We  note  that 
at  the  extreme  temperature  of  6000  °R  (3333.   K)  the  pressure  dependence  is  still  relatively 
small. 

The  uncertainty  in  viscosity  is  estimated  as  follows:    About  2%  for  temperature  up  to 
1080°R,  and  then  gradually  increasing  to  about  5%  at  6000°R. 
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Figure  11.        The  Viscosity  of  Dissociating  Hydrogen  at  Various  Pressures 

7.  7    Thermal  Conductivity 

For  temperatures  up  to  180  °R,the  analytical  representation  of  the  experimental 
results  by  Roder  and  Diller   (1970)  was  chosen  as  source  data  for  para  hydrogen  for  dilute 
gas,    dense  gas,   and  liquid  states.     Dilute  gas  values  for  normal  and  equilibrium  hydrogen 
differ      considerably  from  the  para  values  and  are  calculated  from  the  change  in  the  ideal 
gas  heat  capacity,    C^-.     The  density  dependence  for  n-  and  e-   hydrogen  is  taken  to  be  the 
same  as  para. 

At  higher  temperatures  we  again  select  the  critical  evaluation  by  Hanley,    et  al.,  (1970) 
for  the  dilute  gas  and  the  Enskog  theory  calculation  (Hanley,    et  al.  ,    1972)  for  the  density 
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(pressure)  dependence.     At  180  °R,taking  the  difference  between  normal  and  para  hydrogen 
into  account,    the  agreement  between  the  experimental  results  by  Roder  and  Diller  (1970, 
para),   Golubev  and  Kalsina  (1964,    normal),   and  the  calculation  by  Hanley,    et  al.    (1972)   is 
as  good  as  1%  for     pressures  up  to    1500  psia  and  no  worse  than  4%  at  pressures  of  5000 
psia.     Similar  to  viscosity,  the  dense  gas  contribution  for  thermal  conductivity  vanishes  at  a 
particular  temperature  (1000  K).     At  higher  temperatures  the  potential  function  calculations 
agree  reasonably  well  with  the  data  of  Blais  and  Mann  (I960).     At  yet  higher  temperatures, 
the  effects  of  dissociation  on  thermal  conductivity  become  marked,    see    figure  12.     The 
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Figure  12.      The  Thermal  Conductivity  of  Dissociating  Hydrogen  at  Various  Pressures 


calculations  become  sensitive  to  the  exact  model    or  intermolecular  potential  chosen,    see 
for  example,     figure  9  in  the  paper  of  Grier  (1962).     Two  sets  of  data  exist  for  the  thermal 
conductivity  in  the  region  of  dissociation.     The  theoretical  estimates  selected  (Hanley,    et  al., 
1970)  are  in  excellent  agreement  with  the  shock  wave  experiment  of  Bethringer,    et  al.    (1968). 
The  data  of  Israel,    et  al.,  (1966)  differ     markedly  from  the  former  experiment  and  from  any 
potential  function  model  that  could  be  chosen.     A  possible  explanation  is  that  the   sintered 
material  used  in  this  experiment  enhanced  dissociation.     It  almost  appears  as   if  the 
dissociation  temperatures  are  lowered  approximately  250  K  by  an  apparatus   effect. 
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The  uncertainty  in  thermal  conductivity  is  estimated  to  be  3%  for  temperatures 
between  25  and  180  °R,  increasing  to  6%  near  370  CR;      about  6%  for  temperatures  up  to 
3300  °R  and  then  increasing  gradually  to  10%  at  6000  °R.     Note,   however,    that  thermal 
conductivity  is  expected  to  be  infinite  at  the  critical  point. 

7.  8    Generating  Functions 
The  programs  described  in  the  bulk  of  this  report  are  designed  to  perform  linear 
interpolations  in  arrays  of  data  of  the  various  properties.     Each  property  value  used  in  one 
of  the  arrays  has  to  be  computed  in  a   separate  operation  and  entered  into  the  array.     We  have 
used  several  major  computer  programs  to  accomplish  this  task;  to  distinguish  them  from 
the  linear  interpolation  programs, we  call  the  major  programs  generating  functions.     An 
example  of  one  such  generating  function  is  the  program  VALUES  or  the  subroutine  THERMO 
which  was  developed  in  the  preparation  of  the  tables  presented  in  NBS  Monograph  94  (Roder, 
et  al.  ,    1965).     Since  a  number  of  users  have  requested  copies  of  the  generating  functions, 
we  list  them  below. 

Table  9.     Generating  Functions  * 


Temperature  Range 

Program  Name 

Input 

Output 

25  -   180  °R 

program  VALUES  or 
subroutine  THERMO 

PT 

p,    S,  H,    U,    Cp,  Cv, 

(if)p'  (f?)T« w 

function  CONC 

T,    p,    Cp 

k 

function  VISCOS 

P.    T 

T) 

program  STAT 

T 

o  -   p  composition  and 
ideal  gas  functions, 
thermodynamic 

180°R  and  up 

program  GO 

P,T 

p,    S,H,U,    Cp,    Cv, 
Y,    (J),    k,    n 

*    A  recent  revision  incorporates  the  low  temperature  decks  except  STAT,    an  improved 
equation  of  state  for  high  temperatures,   and  a  new  extrapolation  to  10,000  psia  into 
a  single  program  HYPROP,    see  McCarty  and  Weber  (1972). 
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APPENDIX  A 

Tables  of  Values  for  Selected  Isobars 

Parahydrogen  (Output  from  PT  programs) 38 

Equilibrium  Hydrogen  (Output  from  ET  programs).    .    .  45 
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APPENDIX  B.     PROGRAM  LISTINGS:     CDC  AND  IBM 

CDC  Page  No.  IBM 


Parahydrogen 


PTENTH 

53 

PTENTH 

PTDENS 

57 

PTDENS 

PTENTR 

61 

PTENTR 

PTCOND 

65 

PTCOND 

PTVISC 

68 

PTVISC 

PTSOUN 

72 

PTSOUN 

PTCP,   PTCV, 

PTGAMM 

75 

PTCP,   PTCV, 

PTGAMM 

PTLFAC 

80 

PTLFAC 

TSATH 

83 

TSATH 

PSATH 

84 

PSATH 

PHTEMP 

85 

PHTEMP 

PHDENS 

87 

PHDENS 

PHENTR 

91 

PHENTR 

PHCOND 

94 

PHCOND 

PHVISC 

98 

PHVISC 

PHSOUN 

103 

PHSOUN 

PHCP,   PHCV, 

PHGAMM 

105 

PHCP,   PHCV, 

PHGAMM 

PHLFAC 

111 

PHLFAC 

Equilibrium- 

hydrogen 

ETENTH 

115 

ETENTH 

ETDENS 

119 

ETDENS 

ETENTR 

124 

ETENTR 

ETCOND 

127 

ETCOND 

ETVISC 

131 

ETVISC 

ETSOUN 

134 

ETSOUN 

ETCP,   ETCV, 

ETGAMM 

137 

ETCP,   ETCV, 

ETGAMM 

ETLFAC 
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ETLFAC 

ETSATH 
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ETSATH 

EPSATH 
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EPSA  TH 

EHTEMP 

147 

EHTEMP 

EHDENS 

150 

EHDENS 

EHENTR 

153 

EHENTR 

EHCOND 

157 

EHCOND 

EHVISC 

160 

EHVISC 

EHSOUN 

165 

EHSOUN 

EHCP,   EHCV, 

EHGAMM 

168 

EHCP,   EHCV, 

EHGAMM 

EHLFAC 

173 

EHLFAC 
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APPENDIX  C.     PLOTS  OF  MAXIMUM  INTERPOLATION 
ERROR  AND  ESTIMATED  SOURCE  ERROR 
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